############## 
############## 
############## Models
remove(list = ls())

base::library(conflicted)
base::library(tidyverse)
conflict_prefer("filter","dplyr")
base::library(ggplot2)
base::library(dplyr)
base::library(here)
conflict_prefer("here", "here")
base::library(systemfit)
base::library(stargazer)
base::library(DataCombine)

load(here("Data", "Model_Data.RData"))
glimpse(final_data)

summary(final_data$death_rate)
table(final_data$year)

final_data <- final_data %>% 
  mutate(time_trend = ifelse(year == 2004, 1, NA),
         time_trend = ifelse(year == 2008, 2, time_trend),
         time_trend = ifelse(year == 2012, 3, time_trend),
         time_trend = ifelse(year == 2016, 4, time_trend),
         time_trend = ifelse(year == 2020, 5, time_trend)) 

glimpse(final_data)

table(final_data$state)

test_data <- final_data %>%  filter(!(state == "KY")) %>% 
  filter(!(state == "WV"))

eq1_d <- dem_rep ~ death_rate + unemp_rate + 
  log(defl_pcincome) + 
  share_black + share_hisp + share_asian + share_young + share_elder + 
  log1p(share_college) + log1p(share_manuf) +
  log(tot_pop) + as.factor(rural_urban) + state_partiesdiff +
  state + rep_inc + incumbent_cand + lag_dem_rep + time_trend
eq2_d <- abs_rep ~ death_rate + unemp_rate + 
  log(defl_pcincome) + 
  share_black + share_hisp + share_asian + share_young + share_elder + 
  log1p(share_college) + log1p(share_manuf) +
  log(tot_pop) + as.factor(rural_urban) + state_partiesdiff +
  state + rep_inc + incumbent_cand + lag_abs_rep + time_trend
Syst_death <- list(eq1_d, eq2_d)
model_d <- systemfit(Syst_death, method = "SUR", data = test_data)
summary(model_d)


(estia <- cbind(Estimate = coef(model_d), confint(model_d, level = 0.90)))
length(estia[,1])
estia[2,]
eq1ia <- estia[2, ] #death rates - Democrats
estia[73,]
eq2ia <- estia[73, ] #death rates - Abstention

equations <- c("Democratic/Republican", "Abstention/Republican")
results <- rbind(eq1ia, eq2ia)
results

results <- data.frame(cbind(results, equations))
glimpse(results)

results <- results %>% 
  mutate(equations = factor(equations, 
                            levels = c("Democratic/Republican", 
                                       "Abstention/Republican"))) %>% 
  mutate(Estimate = as.numeric(Estimate)) %>% 
  mutate(X5.. = as.numeric(X5..)) %>% 
  mutate(X95.. = as.numeric(X95..))

fig_sur <- results %>% 
  #filter(variables == "Unemp.Rate" | variables == "Unemp.Rate*Polarization") %>% 
  ggplot(aes()) +
  geom_hline(yintercept = 0, colour = "gray37", lty = 2, linewidth = .75) +
  geom_pointrange(aes(x = equations, y = Estimate, ymin = X5..,
                      ymax = X95..), shape = 15, 
                  lwd = 0.8, position = position_dodge(width = .5)) + 
  theme_bw() + 
  ylab("Log-Ratio Form") + xlab(NULL) + #ylim(-0.1, 0.05) +
  theme(axis.text.x = element_text(hjust = 0.5, size = 13),
        axis.text.y = element_text(size = 10),
        axis.title.y = element_text(size = 13),
        legend.text = element_text(size = 14),
        legend.title = element_text(size = 14),
        legend.position = "bottom",
        legend.direction = "horizontal",
        title = element_text(size = 13)) #+
# ggtitle("A. Coefficient Plot")
fig_sur


